library(readr)
library(rms)
library(bkmr)
library(mice)
library(forestplot)
library(corrplot)
setwd('D:/rCode/论文实战/paper3')
demo <- read_csv("demo.csv")
View(demo)
############################# MI补充数据代码 ###################################
library(tidyverse)
library(gtsummary)
library(haven)
library(skimr)
library(corrplot)
library(bkmr)
library(forestplot)
library(rms)
library(mice)
getwd()

#获取2007至2016年所有研究对象
setwd('D:/rCode/论文实战/nhanes')
demo_e <- read_xpt("2007-2008/Demographics/demo_e.xpt")
demo_f <- read_xpt("2009-2010/Demographics/demo_f.xpt")
demo_g <- read_xpt("2011-2012/Demographics/demo_g.xpt")
demo_h <- read_xpt("2013-2014/Demographics/demo_h.xpt")
demo_i <- read_xpt("2015-2016/Demographics/demo_i.xpt")
demo_all <- dplyr::bind_rows(list(demo_e,demo_f,demo_g,demo_h,demo_i))

#获取所有实验室尿酸数据
Laboratory_e <- read_xpt("2007-2008/Laboratory/BIOPRO_e.xpt")
Laboratory_f <- read_xpt("2009-2010/Laboratory/BIOPRO_f.xpt")
Laboratory_g <- read_xpt("2011-2012/Laboratory/BIOPRO_g.xpt")
Laboratory_h <- read_xpt("2013-2014/Laboratory/BIOPRO_h.xpt")
Laboratory_i <- read_xpt("2015-2016/Laboratory/BIOPRO_i.xpt")
Laboratory_all <- dplyr::bind_rows(list(Laboratory_e,Laboratory_f,Laboratory_g,Laboratory_h,Laboratory_i))

#获取实验室中尿液邻苯二甲酸酯数据
PHTHTE_e <- read_xpt("2007-2008/Laboratory/PHTHTE_e.xpt")
PHTHTE_f <- read_xpt("2009-2010/Laboratory/PHTHTE_f.xpt")
PHTHTE_g <- read_xpt("2011-2012/Laboratory/PHTHTE_g.xpt")
PHTHTE_h <- read_xpt("2013-2014/Laboratory/PHTHTE_h.xpt")
PHTHTE_i <- read_xpt("2015-2016/Laboratory/PHTHTE_i.xpt")
PHTHTE_all <- dplyr::bind_rows(list(PHTHTE_e,PHTHTE_f,PHTHTE_g,PHTHTE_h,PHTHTE_i))

#获取BMI数据
BMI_e <- read_xpt("2007-2008/Examination/bmx_e.xpt")
BMI_f <- read_xpt("2009-2010/Examination/bmx_f.xpt")
BMI_g <- read_xpt("2011-2012/Examination/bmx_g.xpt")
BMI_h <- read_xpt("2013-2014/Examination/bmx_h.xpt")
BMI_i <- read_xpt("2015-2016/Examination/bmx_i.xpt")
BMI_all <- dplyr::bind_rows(list(BMI_e,BMI_f,BMI_g,BMI_h,BMI_i))

#获取血清可替宁数据
cotnal_e <- read_xpt("2007-2008/Laboratory/cotnal_e.xpt")
cotnal_f <- read_xpt("2009-2010/Laboratory/cotnal_f.xpt")
cotnal_g <- read_xpt("2011-2012/Laboratory/cotnal_g.xpt")
cot_h <- read_xpt("2013-2014/Laboratory/cot_h.xpt")
cot_i <- read_xpt("2015-2016/Laboratory/cot_i.xpt")
cot_all <- dplyr::bind_rows(list(cotnal_e,cotnal_f,cotnal_g,cot_h,cot_i)) 

#获取饮酒数据
alq_e <- read_xpt("2007-2008/Questionnaire/alq_e.xpt")
alq_f <- read_xpt("2009-2010/Questionnaire/alq_f.xpt")
alq_g <- read_xpt("2011-2012/Questionnaire/alq_g.xpt")
alq_h <- read_xpt("2013-2014/Questionnaire/alq_h.xpt")
alq_i <- read_xpt("2015-2016/Questionnaire/alq_i.xpt")
alq_all <- dplyr::bind_rows(list(alq_e,alq_f,alq_g,alq_h,alq_i)) 


#获取尿肌酐数据
alb_cr_e <- read_xpt("2007-2008/Laboratory/alb_cr_e.xpt")
alb_cr_f <- read_xpt("2009-2010/Laboratory/alb_cr_f.xpt")
alb_cr_g <- read_xpt("2011-2012/Laboratory/alb_cr_g.xpt")
alb_cr_h <- read_xpt("2013-2014/Laboratory/alb_cr_h.xpt")
alb_cr_i <- read_xpt("2015-2016/Laboratory/alb_cr_i.xpt")
alb_cr_all <- dplyr::bind_rows(list(alb_cr_e,alb_cr_f,alb_cr_g,alb_cr_h,alb_cr_i))

#获取血压问卷数据
bpq.e <- read_xpt("2007-2008/Questionnaire/bpq_e.xpt")
bpq.f <- read_xpt("2009-2010/Questionnaire/bpq_f.xpt")
bpq.g <- read_xpt("2011-2012/Questionnaire/bpq_g.xpt")
bpq.h <- read_xpt("2013-2014/Questionnaire/bpq_h.xpt")
bpq.i <- read_xpt("2015-2016/Questionnaire/bpq_i.xpt")
bpq_all <- dplyr::bind_rows(list(bpq.e,bpq.f,bpq.g,bpq.h,bpq.i))

#获取活动数据
paq_e <- read_xpt("2007-2008/Questionnaire/paq_e.xpt")
paq_f <- read_xpt("2009-2010/Questionnaire/paq_f.xpt")
paq_g <- read_xpt("2011-2012/Questionnaire/paq_g.xpt")
paq_h <- read_xpt("2013-2014/Questionnaire/paq_h.xpt")
paq_i <- read_xpt("2015-2016/Questionnaire/paq_i.xpt")
paq_all <- dplyr::bind_rows(list(paq_e,paq_f,paq_g,paq_h,paq_i))

#获取糖化血红蛋白数据
ghb_e <- read_xpt("2007-2008/Laboratory/ghb_e.xpt")
ghb_f <- read_xpt("2009-2010/Laboratory/ghb_f.xpt")
ghb_g <- read_xpt("2011-2012/Laboratory/ghb_g.xpt")
ghb_h <- read_xpt("2013-2014/Laboratory/ghb_h.xpt")
ghb_i <- read_xpt("2015-2016/Laboratory/ghb_i.xpt")
ghb_all <- dplyr::bind_rows(list(ghb_e,ghb_f,ghb_g,ghb_h,ghb_i))
ghb_all$LBXGH
#获取空腹血糖数据
glu_e <- read_xpt("2007-2008/Laboratory/glu_e.xpt")
glu_f <- read_xpt("2009-2010/Laboratory/glu_f.xpt")
glu_g <- read_xpt("2011-2012/Laboratory/glu_g.xpt")
glu_h <- read_xpt("2013-2014/Laboratory/glu_h.xpt")
glu_i <- read_xpt("2015-2016/Laboratory/glu_i.xpt")
glu_all <- dplyr::bind_rows(list(glu_e,glu_f,glu_g,glu_h,glu_i))
glu_all$LBXGLU
#获取OGTT数据
OGTT_e <- read_xpt("2007-2008/Laboratory/OGTT_e.xpt")
OGTT_f <- read_xpt("2009-2010/Laboratory/OGTT_f.xpt")
OGTT_g <- read_xpt("2011-2012/Laboratory/OGTT_g.xpt")
OGTT_h <- read_xpt("2013-2014/Laboratory/OGTT_h.xpt")
OGTT_i <- read_xpt("2015-2016/Laboratory/OGTT_i.xpt")
OGTT_all <- dplyr::bind_rows(list(OGTT_e,OGTT_f,OGTT_g,OGTT_h,OGTT_i))
OGTT_all$LBXGLT

#获取CKD数据
kiq_u_e <- read_xpt("2007-2008/Questionnaire/kiq_u_e.xpt")
kiq_u_f <- read_xpt("2009-2010/Questionnaire/kiq_u_f.xpt")
kiq_u_g <- read_xpt("2011-2012/Questionnaire/kiq_u_g.xpt")
kiq_u_h <- read_xpt("2013-2014/Questionnaire/kiq_u_h.xpt")
kiq_u_i <- read_xpt("2015-2016/Questionnaire/kiq_u_i.xpt")
kiq_u_all <- dplyr::bind_rows(list(kiq_u_e,kiq_u_f,kiq_u_g,kiq_u_h,kiq_u_i))
kiq_u_all$KIQ022
output <- plyr::join_all(list(demo_all, Laboratory_all,PHTHTE_all,BMI_all,cot_all,alq_all,alb_cr_all,bpq_all
                              ,paq_all,ghb_all,glu_all,OGTT_all,kiq_u_all),
                         by='SEQN',type='full')
#年龄筛选
output <- output[output$RIDAGEYR>=20,]
dim(output)
#尿酸数据筛选
output <- output[!is.na(output$LBDSUASI),]
dim(output)
#邻苯二甲酸酯数据筛选
output <- output[!is.na(output$URXCNP),]
output <- output[!is.na(output$URXCOP),]
output <- output[!is.na(output$URXECP),]
output <- output[!is.na(output$URXMBP),]
output <- output[!is.na(output$URXMC1),]
output <- output[!is.na(output$URXMEP),]
output <- output[!is.na(output$URXMHH),]
output <- output[!is.na(output$URXMOH),]
output <- output[!is.na(output$URXMZP),]
output <- output[!is.na(output$URXMIB),]
#BMI数据筛选
output <- output[!is.na(output$BMXBMI),]#8409
dim(output)
#家庭贫困数据筛选
output <- output[!is.na(output$INDFMPIR),]#7650
dim(output)
#饮酒数据筛选 
output <- output[!is.na(output$ALQ101),]
outnadrink1 <- which(is.na(output$ALQ101))
dim(output)
#血清可替宁数据筛选
output <- output[!is.na(output$LBXCOT),]#8407
#教育程度数据筛选
output <- output[!is.na(output$DMDEDUC2),]#7047
#尿肌酐数据筛选
output <- output[!is.na(output$URXCRS),]#7045
dim(output)
#高血压数据筛选
a <- output
output <- output[!is.na(output$BPQ020),]#7045
dim(a)
#体力活动数据筛选PAQ665
#output$sports <- ifelse(output$PAQ650==2,0,ifelse(output$PAQ665==2,1,NA))
#output <- output[!is.na(output$PAQ665),]#7045
#dim(output)
#判断糖尿病
demo1 <- output
demo1$diabetes <- ifelse(demo1$LBXGH>5.7 | demo1$LBXGLU>100 | demo1$LBXGLU>140,
                         1,0)
#性别
demo1$Sex <- ifelse(demo1$RIAGENDR==1,'male','female')
#bmi
demo1$BMI <- ifelse(demo1$BMXBMI<=24.9,'≤ 24.9', '> 24.9')
# 种族
demo1$race <- recode_factor(demo1$RIDRETH1, 
                            `1` = 'Mexican American',
                            `2` = 'Other Hispanic',
                            `3` = 'White',
                            `4` = 'Black',
                            `5` = 'Other Race')
# 教育程度
demo1$edu <- recode_factor(demo1$DMDEDUC2, 
                           `1` = 'Less than high school',
                           `2` = 'Less than high school',
                           `3` = 'High school graduate/GED or equivalent',
                           `4` = 'College or above',
                           `5` = 'College or above')
# 贫困程度
demo1$FPL <- ifelse(demo1$INDFMPIR<1.0,'< 1.0',
                    ifelse((demo1$INDFMPIR>=1.0&demo1$INDFMPIR<2.0) , "1.0–2.0" ,
                           ifelse((demo1$INDFMPIR>=2.0&demo1$INDFMPIR<4.0) , "2.0–4.0" ,'> 4.0')))
#活动
demo1$activity <- ifelse(demo1$PAQ605==1,'vigorous',
                         ifelse((demo1$PAQ620==1) , "moderate" ,'no'))
# 吸烟 1 是 
demo1$smoke <- ifelse(demo1$LBXCOT>10,'yes','no')
# 饮酒
demo1$drink <- ifelse(demo1$ALQ101==1,'yes','no')
#高血压
demo1$Hypertension <- ifelse(demo1$BPQ020==1,'yes','no')
# ckd 
demo1$CKD <- ifelse(demo1$KIQ022==1,'yes','no')
# 尿肌酐
demo1$urinecreatinine <- demo1$URXCRS
# 年龄
demo1$Age <- demo1$RIDAGEYR
demo1$Hyperuricemia <- ifelse((demo1$Sex=='male' & demo1$LBDSUASI>=416.0),'Yes',
                              ifelse((demo1$Sex=='male' & demo1$LBDSUASI<416.0),'No',
                                     ifelse((demo1$Sex=='female' & demo1$LBDSUASI>=357.0),'Yes',
                                            ifelse((demo1$Sex=='female' & demo1$LBDSUASI<357.0),'No',""))))
demo <-demo1[,c('SEQN', 'RIDAGEYR', 'Sex','Age', 'BMI','race','edu','FPL', 'activity','smoke','drink','Hypertension','diabetes','CKD',
                'urinecreatinine','URXECP','URXMBP','URXMHH','URXMOH','URXMIB','URXCNP','URXCOP','URXMC1','URXMEP','URXMZP','LBDSUASI','Hyperuricemia','URXCRS')]
demoMI <- mice(demo, #数据集
               method = "rf", #采用随机森林插补
               m=5, # 5次插补
               printFlag = FALSE #不显示历史记录
)
demoMI <- complete(demoMI, action = 2)
colnames(demoMI)
demoMI$Hyperuricemiaint <- ifelse(demoMI$Hyperuricemia=='Yes',1,0)
demoMI <- demoMI[,c('Age','Sex','race','edu','FPL','activity','smoke',
                    'drink','Hypertension','CKD','urinecreatinine','BMI','URXECP', 'URXMBP', 
                    'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP',
                    'Hyperuricemia','Hyperuricemiaint','LBDSUASI','URXCRS')]
colnames(demoMI) <- c('Age','Sex','race','edu','FPL','activity','smoke',
                      'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP','MnBP','MEHHP',
                      'MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP',
                      'Hyperuricemia','Hyperuricemiaint','LBDSUASI','URXCRS')
demo <- demoMI
setwd('D:/rCode/论文实战/paper3')
write.csv(demoMI,'demo.csv')







################################## table S3 ####################################
#参与的变量
myVars1 <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
nonvar1 <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
tab1 <- CreateTableOne(vars = myVars1, strata =c("Hyperuricemia")  , data = demo)
tabMat1 <- print(tab1, nonnormal = nonvar1,test = TRUE)
tabMat1 <- tabMat1[,c(-4)]
tabMat1
demo$MiNP


######################## 邻苯二甲酸酯相关性热图 图2 ############################
demo2 <- demo[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')]
data<-cor(demo2,method="spearman")

corrplot(data, method = "number", type = "upper",order = 'hclust',
         tl.col = "black", tl.cex = 0.8, tl.srt = 45,tl.pos = "lt")

corrplot(data, method = "circle", type = "lower",order = 'hclust',
         tl.col = "n", tl.cex = 0.8, tl.pos = "n",
         add = T)



################################## RCS 图S2 ####################################
demoS2 <- demo
demoS2[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')] <- 
  log(demoS2[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')])

ddist<-datadist(demoS2)
options(datadist="ddist")
RCSX <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
RCSNK <- c(3,3,3,3,3,5,3,3,3,3)
i <- 1
par(mfrow=c(2,5))
for (i in 1:10) {
  RCSX1 <- RCSX[c(-i)]

  paste0(RCSX1,collapse = '+')
  fml <- as.formula(paste0('Hyperuricemiaint ~rcs(',RCSX[c(i)],',nk=',RCSNK[c(i)],')+',paste0(RCSX1,collapse = '+'),
                '+Age + Sex + race + edu + FPL + activity + smoke + drink + Hypertension + CKD + urinecreatinine + BMI'))
  fit <-lrm(fml,data=demoS2,x=FALSE)
  pred_OR<-Predict(fit,name=RCSX[c(i)],ref.zero=TRUE,fun=exp)

  ylim.bot<-min(pred_OR[,"lower"])
  ylim.top<-max(pred_OR[,"upper"])
    plot(pred_OR[,1],pred_OR[,"yhat"], 
         xlab = RCSX[c(i)],ylab = paste0("OR"),
         type = "l",ylim = c(ylim.bot,ylim.top),
         col="red",lwd=2)+
      lines(pred_OR[,1],pred_OR[,"lower"],lty=2,lwd=1.5)+
      lines(pred_OR[,1],pred_OR[,"upper"],lty=2,lwd=1.5)+
      lines(x=range(pred_OR[,1]),y=c(1,1),lty=3,col="grey40",lwd=1.3) +
      points(1,1,pch=16,cex=1.2)#需要设置参考值,
}



################################## RCS 图S3 ####################################
demoS2 <- demo
demoS2[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')] <- 
  log(demoS2[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')])

ddist<-datadist(demoS2)
options(datadist="ddist")
RCSX <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
RCSNK <- c(3,3,3,3,3,5,3,3,3,3)
i <- 1
par(mfrow=c(2,5))
for (i in 1:10) {
  RCSX1 <- RCSX[c(-i)]
  
  paste0(RCSX1,collapse = '+')
  fml <- as.formula(paste0('LBDSUASI ~rcs(',RCSX[c(i)],',nk=',RCSNK[c(i)],')+',paste0(RCSX1,collapse = '+'),
                           '+Age + Sex + race + edu + FPL + activity + smoke + drink + Hypertension + CKD + urinecreatinine + BMI'))
  fit <-lrm(fml,data=demoS2,x=FALSE)
  pred_OR<-Predict(fit,name=RCSX[c(i)],ref.zero=TRUE,fun=exp)
  
  ylim.bot<-min(pred_OR[,"lower"])
  ylim.top<-max(pred_OR[,"upper"])
  plot(pred_OR[,1],pred_OR[,"yhat"], 
       xlab = RCSX[c(i)],ylab = paste0("OR"),
       type = "l",ylim = c(ylim.bot,ylim.top),
       col="red",lwd=2)+
    lines(pred_OR[,1],pred_OR[,"lower"],lty=2,lwd=1.5)+
    lines(pred_OR[,1],pred_OR[,"upper"],lty=2,lwd=1.5)+
    lines(x=range(pred_OR[,1]),y=c(1,1),lty=3,col="grey40",lwd=1.3) +
    points(1,1,pch=16,cex=1.2)#需要设置参考值,
}



############################## bkmr单变量暴漏 图S4 #############################
demo <- read_csv("demo.csv")
demoS5 <-demo[,c('Age','Sex','race','edu','FPL','activity','smoke',
                 'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP',
                 'MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP',
                 'Hyperuricemiaint','LBDSUASI')]
demoS5 <- demoS5[c(1:200),]
y <- demoS5$Hyperuricemiaint
Z <- demoS5[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP',
               'MCOP','MCPP','MEP','MBzP')]
X <- demoS5[,c('Age','Sex','race','edu','FPL','activity','smoke',
               'drink','Hypertension','CKD','urinecreatinine','BMI')]
for(i in 1:10){
  Z[,i] <- log(Z[,i])
}
set.seed(2012)
fitkm <- kmbayes(y, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'binomial',est.h = TRUE)
pred.resp.univar <- PredictorResponseUnivar(fit = fitkm)
ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity") + 
  facet_wrap(~ variable) +
  ylab("h(expos)") +
  xlab("Phthalate metabolites") +
  xlim(-2,6)



###################### bkmr三个分位数的暴露-应答趋势 图S5 ######################
demo <- read_csv("demo.csv")
demoS4 <-demo[,c('Age','Sex','race','edu','FPL','activity','smoke',
                 'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP',
                 'MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP',
                 'Hyperuricemiaint','LBDSUASI')]
demoS4 <- demoS4[c(1:200),]
y <- demoS4$Hyperuricemiaint
Z <- demoS4[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP',
               'MCOP','MCPP','MEP','MBzP')]
X <- demoS4[,c('Age','Sex','race','edu','FPL','activity','smoke',
               'drink','Hypertension','CKD','urinecreatinine','BMI')]
for(i in 1:10){
  Z[,i] <- log(Z[,i])
}
set.seed(2012)
fitkm <- kmbayes(y =y, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'binomial',est.h = TRUE)
pred.resp.bivar <- PredictorResponseBivar(fit =fitkm , min.plot.dist = 1)
pred.resp.bivar.levels <- PredictorResponseBivarLevels(  
  pred.resp.df = pred.resp.bivar, 
  Z = Z, qs = c(0.1, 0.5, 0.9))
ggplot(pred.resp.bivar.levels, aes(z1, est)) + 
  geom_smooth(aes(col = quantile), stat = "identity") + 
  facet_grid(variable2 ~ variable1) +
  ggtitle("h(expos1 | quantiles of expos2)") +
  xlab("expos1")



##################### bkmr三个分位数的单次暴露风险比较 图S6 ####################
demo <- read_csv("demo.csv")
demoS4 <-demo[,c('Age','Sex','race','edu','FPL','activity','smoke',
                 'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP',
                 'MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP',
                 'Hyperuricemiaint','LBDSUASI')]
demoS4 <- demoS4[c(1:200),]
y <- demoS4$Hyperuricemiaint
Z <- demoS4[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP',
               'MCOP','MCPP','MEP','MBzP')]
X <- demoS4[,c('Age','Sex','race','edu','FPL','activity','smoke',
               'drink','Hypertension','CKD','urinecreatinine','BMI')]
for(i in 1:10){
  Z[,i] <- log(Z[,i])
}
set.seed(2012)
fitkm <- kmbayes(y =y, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'binomial',est.h = TRUE)
risks.singvar <- SingVarRiskSummaries(fit = fitkm, y = y, Z = Z, X = NULL, 
                                      qs.diff = c(0.25, 0.75),
                                      method = "exact")
ggplot(risks.singvar, aes(variable, est, ymin = est - 1.96*sd, 
                          ymax = est + 1.96*sd, col = q.fixed)) + 
  geom_pointrange(position = position_dodge(width = 0.75)) + 
  coord_flip()



############################# bkmr总体效应 图S7B ##############################
demo <- read_csv("demo.csv")
demoS4 <-demo[,c('Age','Sex','race','edu','FPL','activity','smoke',
                'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP',
                'MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP',
                'Hyperuricemiaint','LBDSUASI')]
demoS4 <- demoS4[c(1:200),]
y <- demoS4$Hyperuricemiaint
Z <- demoS4[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP',
              'MCOP','MCPP','MEP','MBzP')]
X <- demoS4[,c('Age','Sex','race','edu','FPL','activity','smoke',
              'drink','Hypertension','CKD','urinecreatinine','BMI')]
for(i in 1:10){
  Z[,i] <- log(Z[,i])
}
set.seed(2012)
fitkm <- kmbayes(y =y, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'binomial',est.h = TRUE)
risks.overall = OverallRiskSummaries(fit=fitkm,qs=seq(0.25,0.75,by=0.05),q.fixed = 0.5)
risks.overall
ggplot(risks.overall,aes(quantile,est,ymin=est-1.96*sd,ymax=est+1.96*sd))+
  geom_hline(yintercept = 0,lty=2,col='brown')+
  geom_pointrange()



################################# 森林图 图3 ###################################
demo3 <- demo
demo3$MECPP_quantile <- cut(demo3$MECPP,breaks = quantile(demo3$MECPP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))
demo3$MnBP_quantile <- cut(demo3$MnBP,breaks = quantile(demo3$MnBP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MEHHP_quantile <- cut(demo3$MEHHP,breaks = quantile(demo3$MEHHP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MEOHP_quantile <- cut(demo3$MEOHP,breaks = quantile(demo3$MEOHP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MiBP_quantile <- cut(demo3$MiBP,breaks = quantile(demo3$MiBP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$cxMiNP_quantile <- cut(demo3$MiNP,breaks = quantile(demo3$MiNP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MCOP_quantile <- cut(demo3$MCOP,breaks = quantile(demo3$MCOP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MCPP_quantile <- cut(demo3$MCPP,breaks = quantile(demo3$MCPP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MEP_quantile <- cut(demo3$MEP,breaks = quantile(demo3$MEP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MBzP_quantile <- cut(demo3$MBzP,breaks = quantile(demo3$MBzP),labels = c('Q1', 'Q2', 'Q3', 'Q4')) 
varis<- c('MECPP_quantile','MnBP_quantile','MEHHP_quantile',
          'MEOHP_quantile','MiBP_quantile','cxMiNP_quantile','MCOP_quantile',
          'MCPP_quantile','MEP_quantile','MBzP_quantile')
sinLogistic <- function(x){
  my_list<-list(c(1))
  for (variable in x) {
    len <- length(my_list)
    f<- as.formula(paste0('Hyperuricemiaint ~ ',variable,'+urinecreatinine'))
    glmMECPP<-glm(f,data = demo3,family =  binomial )
    my_list[[len+1]] <- glmMECPP
  }
  return(my_list)
}
allSin <-sinLogistic(varis)

my_listF <- list(c(1))
for (variable in c(2:11)) {
  print(variable)
  len <- length(my_listF)
  fit.result<-summary(allSin[[variable]])
  df1<-fit.result$coefficients
  df2<-confint(allSin[[variable]])
  df3<-cbind(df1,df2)
  df4<-data.frame(df3[-1,c(1,4,5,6)])
  df4$Var<-rownames(df4)
  colnames(df4)<-c("OR","Pvalue","OR_1","OR_2","Var")
  df5<-df4[,c(5,1,2,3,4)]
  df5$OR_mean<-df5$OR
  df5$OR<-paste0(round(df5$OR,2),
                 "(",
                 round(df5$OR_1,2),
                 "~",
                 round(df5$OR_2,2),
                 ")")
  df5$Pvalue<-round(df5$Pvalue,3)
  my_listF[[len+1]] <- df5
}
allfros <- dplyr::bind_rows(list(my_listF[[2]],my_listF[[3]],my_listF[[4]],my_listF[[5]],my_listF[[6]],my_listF[[7]],my_listF[[8]],my_listF[[9]],my_listF[[10]],my_listF[[11]]))
fp<-allfros
forestplot(labeltext=as.matrix(fp[,1:3]),
           mean=fp$OR_mean,
           lower=fp$OR_1,
           upper=fp$OR_2,
           zero=0.25,
           boxsize=0.2,
           graph.pos=2)












